{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "streaming-filename",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical integration using quad()\n",
    "*March 19, 2021*\n",
    "\n",
    "- First import some required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "checked-marker",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import quad"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plastic-coalition",
   "metadata": {},
   "source": [
    "An integral that comes up when trying to find the chemical potential of a system of identical Bosons is:\n",
    "\n",
    "$\\int_0^\\infty \\dfrac{x^2}{e^{x^2-\\eta}-1}\\, dx$\n",
    "\n",
    "where $\\eta=\\mu/k_\\mathrm{B}T$ and $\\mu$ is the chemical potential.\n",
    "\n",
    "This short tutorial deminstrate how to use *quad()* to numberically evaluate this integral for the $\\eta = 0$ case.\n",
    "\n",
    "\n",
    "- First, define a function that represents the integrand of the integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "floppy-inspection",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrand(x):\n",
    "    return x**2/(np.exp(x**2) - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "damaged-montgomery",
   "metadata": {},
   "source": [
    "- Next, we call the function inside *quad()* from the SciPy module to perfrom the numerical integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "verbal-movie",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-27-472f14518a59>:2: RuntimeWarning: overflow encountered in exp\n",
      "  return x**2/(np.exp(x**2) - 1)\n"
     ]
    }
   ],
   "source": [
    "y, err = quad(integrand, 0, np.inf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "superb-amino",
   "metadata": {},
   "source": [
    "- *quad()* returns two values.  The first (y) is the value of the integral and the second (err) is an estimate of the error in the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "original-recorder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1575786866969464, 8.140518531748718e-09)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y, err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "representative-welcome",
   "metadata": {},
   "source": [
    "- Notice also that we got and overflow warning.  This occurs because when $x$ is very large $e^{x^2}$ becomes exceedingly large.  If we want, we can avoid this warning by simply truncating the integration using some maximum value of $x$, say $\\infty\\to 10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "breathing-crown",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.1575786866970472, 7.989498536842354e-12)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y, err = quad(integrand, 0, 10)\n",
    "y, err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "operating-rugby",
   "metadata": {},
   "source": [
    "For this choise of the upper integration limit, there was no warning and our result for y changed by only a completely negligible amount."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}